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ABSTRACT 

We search for the signature of an isotropic stochastic gravitational-wave background 
in pulsar timing observations using a frequency-domain correlation technique. These 
observations, which span roughly 12 yr, were obtained with the 64-m Parkes radio 
telescope augmented by public domain observations from the Arecibo Observatory. A 
wide range of signal processing issues unique to pulsar timing and not previously pre- 
sented in the literature are discussed. These include the effects of quadratic removal, 
irregular sampling, and variable errors which exacerbate the spectral leakage inher- 
ent in estimating the steep red spectrum of the gravitational- wave background. These 
observations are found to be consistent with the null hypothesis, that no gravitational- 
wave background is present, with 76 percent confidence. We show that the detection 
statistic is dominated by the contributions of only a few pulsars because of the inho- 
mogeneity of this data set. The issues of detecting the signature of a gravitational- wave 
background with future observations are discussed. 
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1 INTRODUCTION 

Times-of-arrival (ToAs) of high signal-to- noise (S/N) ratio 
integrated pulses from millisecond pulsars (MSPs) can be 
measured very precisely, often with sub-/is uncertainties. 
The rotational stability of a MSP implies that a simple 
model of the pulsar can be developed to make accurate pre- 
dictions of these ToAs. Comparing the measured ToAs with 
these predictions enables the study of many astrophysical 
phenomena; for example, this process le d to evidence for the 
existence of gravitational waves (GWs; iTavlor fc Weisberd 
Ll982). As more MSPs are discovered and instrumentation 
is improved, it is becoming likely that pulsar observations 
will lead to the direct detection of G Ws, using the ir ef- 
fect on ToAs described independently bv lSazhinI l|l978l ) and 

* E-mail: dyardley@physics.usyd.edu.au (DRBY) 



iDetweileJ ((1971). A GW will cause a perturbation in the 
ToA when it passes the pulsar and again when it passes 
the Earth. The perturbations that would be detectable with 
pulsar timing are expected to have amplitu des of ~10ns 
and t imescales greater than one year (see, e.g.. lSesana et af] 
I2OO9I) . 

The P arkes Pulsar Timing Array (PPTA) project (e.g., 
IVerbiest et al. 2010) is using the Parkes radio telescope and 
advanced instrumentation to time 20 MSPs over a period 
of at least 5 years. With careful calibration and long inte- 
grations, the majority of the pulsars are yielding weighted 
root-mean-square (rms) residuals below 1 us, with a few be- 
low 200ns ((Manchester 2010). While some of the PPTA 
pulsars do show "timing noise" (low- frequency timing in- 
st abilities which are unmodelled by conventional analyses), 
IVerbiest et al.l (|2009l ') showed that this will not prohibit GW 
detection with the PPTA pulsars. The project will allow ex- 
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amination of correlated signals between the different pulsars, 
includ ing detecting v a riatio n s in the terrestrial timescale 
(e.g., IPetit fc TaveUal 1 19961 : iHobbs et aP '2010b'), detect- 
ing e rrors in the Solar-System ephemeris ( Champ ion et al.l 
l201(]f ). and providing constraining limits on, or a detection 
of, low-frequency GWs. The project has been ongoing since 
late 2004. Observations of some of the PPTA pulsars have 
been made at the Parkes observatory since 1994, albeit with 

less regularity and pre cision. 

Recent work jYardley et al.l I2OI0I: 
2OIOI: Ivan Haasteren fc LevinI 120101: 



Burt et al.l 

n=i ISesana et all |2009|: 

Corbin fc Cornistil 2010l ) has addressed the detectability 



of individual sources of GWs in pulsar timing residuals 
and shows that it is unlikely that current instrumentation 
will allow a detection. However, if the universe contains 
many such sources of GWs, these sources will form an 
isotrop i c stoc hastic gravitational- wave background (GWB). 
ISazhinI ll 19781). ipetweilei] (|1979( ). iHellings fc DownsI (| 19831 ) 
and Jenet et al.l ( 20051 ) have described how pulsar timing 
arrays (PTAs), such as the PPTA, can directly detect 
such a background of ~nHz frequency GWs. For each 
pulsar, this GWB would cause ToA perturbations that are 
correlated between pulsar pairs in a quadrupolar fashion. 
This correlation, which depends only on the angle between 
the p air of pulsars as shown in Figure 1 (|Hellings fc Downg 
I1983D . provides an unambiguous signature of the GWB. 
The functional form of this signature is given by: 



C(f»j) = ^x\ogx - I + ^ 



(1) 



where a; = [1 — cos{9ij)]/2 and 6ij is the angle between pul- 
sars i and ,7 subtended at the observer (jHellings fc Downg 
ll983l : IJenet et akllioosl ). The function (^{Oij) is independent 
of GW frequency, and is derived assuming GWs are de- 
sc ribed by g e neral relativity; other GW modes are analysed 
m iLee et all l|2008l ) but are not considered in this paper. We 
believe that a first detection of the GWB is only possible 
via an unambiguous detection of this expected correlation. 
In view of the widespread interest in such a detection, we 
have designed a detection procedure that can show this sig- 
nature in an easily discernible and convincing manner. 

Several techniques have alre ady been propose d in the 
literature to both limit |Romani fc Taylor 1983: Kaspi et al 



19941 : iThorsett fc Dewevlll996l: lLommenll2002l: Ijenet et al 



20od ) and detect (| Jenet et al I BoOSl : I Anholm et al.l |2009| ) 
the GWB. However, these methods have not taken into ac- 
count all the details of optimally treating pulsar timing data, 
or are restricted to pa rticular observations. Appl ication of 
a Bayesian technique (|van Haasteren et al.l |2009| ) is ongo- 
ing work, and our method is completely independent. The 
lowest published li mit for a GWB cau sed by supermassive 
binary black holes (|jenet et al.l I2OO6I ) beg ins to constrain 
the p arameters of galaxy evolution (e.g., Wvithe fc Loebl 
l2003l ). cosmic strings (e.g., iDamour fc Vilenkinll2005l ) and 



relic GWs from the Big Bang (e.g., Maggior j 2000l ) . Further 
improvements in sensitivity could either enable detection of 
GWs or rule out most proposed models of a GWB. 

The GWB de tection techn i que w e present here is based 
on the method of I Jenet et all l|2005l) . It improves on their 
technique in a number of ways: 

- we study the pairwise correlation described by 
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Figure 1. The expected correlation in pulsar timing residuals 
due to an isotropic stochastic GWB. The abscissa gives the angle 
subtended at the observer by a particular pulsar pair. The ordi- 
nate gives the expected correlation between the timing residuals 
of that pair. This signal is independent of the GW frequency and 
assumes that GWs behave as predicted by general relativity. 



iHellings fc DownsI l|l983l l in the form of pairwise cross power 
spectra; 

- we obtain independent estimates of the GWB from each 
frequency component in each cross power spectrum; 

- we use an optimally weighted linear combination of the 
cross power estimates as the detection statistic; 

- we account for the effect of different overlapping times- 
pans between the pulsar pairs; 

- we calibrate the cross power spectra and their estimated 
errors using simulations that completely account for the fit- 
ting of the pulsar timing model. 

Our technique is not optimal for bounding the GWB 
with these observations because the variation in S/N ra- 
tio between pulsars is too large. This means that there are 
not enough significant cross power spectra to compensate 
for the low value of the average cross correlation. A tighter 
bound on the GWB amplitude could be obtained with these 
observations using the amplitudes o f the individual power 
spectra (similar to I Jenet et al.|[2006l l. However, a detection 
algorithm cannot be based on the amplitudes of individual 
power spectra because there are many unknown contribu- 
tions to those power spectra. We discuss a n u mber of is- 
sues that are common to both the lJenet et al.l (|2006l ) limit 
technique and any limit technique based on measuring the 
GWB-induced correlation between pulsars. Such issues in- 
clude the estimation of power spectra when the sampling 
is irregular and the ToA uncertainties are variable, and the 
effects of fitting the timing model. 

In !j2]we describe the observations and the analysis that 
led to the timing residuals we use in this paper. J^jdescribes 
the theoretical background and our method for making a 
detection of the isotropic stochastic GWB. [|4] describes the 
results obtained, !jS]describes their implications and the out- 
standing issues for GWB detection via pulsar timing, and ^6] 
summarises our conclusions. 
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Table 1. Basic information for the Verbiest et al. (2008, 2009) data sets. 
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^ There is a gap of --^ 11 years between the end of the observations presented in lKaspi et al. 1 1I1994I ') and the beginning of observations with 
the Parkes telescope. In our analysis we use the Arecibo observations of PSR J1857+0943 only to assist in the estimation of the pulsar 
parameters and then discard the Arecibo residuals in further processing. 

We have altered the v alue of the phase offsets between different observing systems for these timing residuals compared with the analysis 

of lVerbiest et al ] l|2009l 'l. which lowers our measured rms. 

'-■ This time series features several large gaps and incl udes thelKaspi et al" I lll994l'l data. 

* These values differ slightly from those presented in i Verbiest et al.l l|2009l )"because we have removed duplicated observations in five pulsars, 
and corrected a minor processing error involving the uncertainties on observations made with different observing systems. 



2 OBSERVATIONS 

The 20 pulsars used in this paper were observed for ~10 min 
to 1 h in each observation, depending on the hardware being 
used at the time. Since 2005, the typical integration time on 
most pulsars is ~ 1 h. For each observation a mean pulse pro- 
file was formed using an ephemeris which "folds" the data 
at the apparent pulse period. Observations of each pulsar 
were made every few weeks (although there are some gaps 
of many months) and the observations span many years, as 
shown in Table [T] and Figure [2] The time shift between a 
standard pulse profile and the ob served profile is measured 
using the technique described in iTavloil (|l992l ). as imple- 
mented within the PAT rout i ne of the PSRCHIVE so ftware 
package (jHotan et al.ll2004l : Ivan Straten et al.ll2010l ). This 
measurement results in an estimate of a ToA and its un- 
certainty. The observatory timescale was referenced to Uni- 
versal Coordinated Time and post-corrected to Terrestrial 
Time as realised by International Atomic Time, abbrevi- 
ated to TT(TAI). The effect of corrections to this timescale 
published by the Bureau International des Poids et Mesures 
(BIPM) is discussed in more detail in i ]5.3l The ToAs were 
transformed to a barycen tric arrival tim e using the DE405 
Solar-System ephemeris (|Standishl |2004| ). The barycentric 
ToAs were then fit ted with a timing model using the tempo 2 
software package (|Hobbs et al ] |2006l : lEdwards et alll2006l ). 
We refer to the differences between the observed and pre- 
dicted ToAs as the "timing residuals". Statistically signif- 
icant timing residuals represent physical effects that have 



not been included in the timing model, and can have many 
causes, such as incomplete polarisation calibration, timing 
noise intrinsic to the pulsar system, fluctuations in the ISM, 
errors in the Solar-System ephemeris, errors in the terres- 
trial timescale and GWs. More on the techniques of pul- 
sar timing can be fo und in iLorimer fc Krameil (|2004l ) and 
lEdwards etal] (120061 ). 

In this paper we u se th e timing resid u als p resented 
bv IVerbiest et al.l l|2008l ) and IVerbiest et al.l (|2009l ). These 
residuals are assembled from a large number of different ob- 
servations with different receivers and even different obser- 
vatories. The observations come primarily from the Parkes 
radio telescope and most were made as part of the PPTA 
project. They are augmented by earlier Parkes observations 
and publicly available observations of PSR J1857-I-0943 and 
PSR J1939+2134 t aken with the Arecibo radio telescope and 
described in iKaspi et al.i (,1994 ). The Arecibo observations of 
PSR J1857-I-0943 were carried out at ~1400MHz and span 
seven years. The Arecibo observations of PSR J1939-I-2134 
were carried out at ~ 1400 MHz and ~2400MHz and span 
eight years. 

The Verbiest et al. (2008, 2009) observations were per- 
formed in the 20 cm (1400 MHz) band, except for PSR 
J0613— 0200 for which a better timing solution was obtained 
in the 50 cm (685 MHz) band. The observations have not 
been fully corrected for variations in the pulse dispersion 
measure (DM). Observations in the 20cm band between 
1994 and November 2002 were made with either one or 
two 128 MHz-wide bands, but these data vary greatly in 
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Figure 2. The timing residuals for the four most influential 
pulsars used in this analysis. The length of the vertical line on 
the left-hand side indicates 10 /xs. The right-hand column gives 
the pulsar's J-name. Noise levels vary significantly both between 
pulsars and at different epochs. 



quality. Observations after November 2002 were taken over 
two 64 MHz- wide observing bands centred at 1341 MHz and 
1405 MHz. For full details of ToA m easurement and data 
processing, see IVerbiest et aD ()2009l ). As mentioned in the 
footnote to Table [1] we have made minor corrections to the 
Verbiest et al. (2008, 2009) observations, though we have not 
performed the full data reduction process already described 
and performed in Verbiest et al. (2008, 2009). We have also 
altered the value of the phase offsets between different ob- 
serving systems for the PSR J1939-I-2134 timing residuals, 
which reduces the measured rms. A summary of the data sets 
is given in Table[T] Of the 20 time series, the four most influ- 
ential are plotted in Figure (2] the ti ming residuals from a ll 
observations are shown in Figure 1 of IVerbiest et al.l (|2009l ). 

The observations were made with a number of different 
observing systems — both the frontend receivers and the 
backend instrumentation have varied over time. Arbitrary 
phase offsets have been fitted for and removed between the 
To As from each different system for a given pulsar. This re- 
duces the noise level in the timing residuals for that pulsar, 
especially over long timescales. These residuals also have a 
number of features which complicate the time series anal- 
ysis and spectral estimation. While the timing residuals of 
most of the pulsars are "white" (i.e., their power spectra are 
independent of frequency) , nine out of the twenty pulsars ex- 
hibit non- white noise. This was determined using a simple 
two-point correlation analysis to determine the degree of cor- 
relation between adjacent residuals using the CHECKWHITE 
plugin to TEMP02. 

The data spans vary widely, ranging from 2.8 years for 
PSR J1824-2452 to 23.3 years for PSR J1939+2134. The 
weighted rms residual also varies over two orders of mag- 
nitude, from 170 ns for PSR J1909-3744 to 15 fis for PSR 
J1939-I-2134. The residuals are also sampled irregularly and 
the sampling is different between pulsars. The ToA uncer- 
tainties for a given pulsar vary widely over short and long 
timescales. This is normally caused by scintillation in the in- 
terstellar medium and upgrades in the receiver and backend 
systems, respectively. In some cases, the magnitude of the 
ToA error bar changes discontinuously during the time series 
due to these upgrades in the observing hardware at Parkes. 
The upgrade which had the largest effect on the quality of 



Table 2. Pulsars with non-stationary timing residuals. For these 
pulsars, we estimate the unweighted rms of the residuals before 
and after an important hardware change at the telescope. 



PSRJ 


Type of 
change 


Epoch 
(MJD) 


RMS before 
change (fis) 


RMS after 
change {fis) 


J1600-3053 


backend 


52654.0 


9.61 


1.31 


J1713-I-0747 


backend 


52462.5 


1.24 


0.48 


J1732-5049 


backend 


52967.5 


7.57 


4.03 


J1744-1134 


backend 


52462.6 


1.54 


1.29 


J2124-3358 


backend 


52984.5 


9.74 


4.64 


J2129-5721 


receiver 


51410.0 


5.47 


3.48 


J2145-0750 


backend 


52975.5 


4.14 


3.17 



the timing residuals was the transition from the Caltech 
incoheren t autocorrelati on spectrometer fast pulsar timing 
machine (|Navarrol 19941) to the Caltech-Parkes-Swinburne 



Recorder 2 ( Bailee 



20031 ). a coherent dedispersion system. 



in late 2002. We therefore attempt to reduce the huge vari- 
ation in the magnitude of the ToA uncertainties so that, in 
subsequent weighted estimates using the timing residuals, 
the weights are spread more evenly across the data set. We 
provide in Table [5] a list of the pulsars for which we have 
calculated the sample variance of the residuals in two dif- 
ferent sections of the time series because of a step-change in 
the quality of the timing residuals. These sample variances 
are added in quadrature with the original error bars in each 
portion before commencing any further processing. For all 
other pulsars there was no significant change in data qual- 
ity at the epoch of the hardware change. We thus calculate 
the sample variance of the whole time series and add it in 
quadrature with the original error bars before any further 
processing. 



3 METHOD 

For all pulsars, the GWB will induce timing residuals with 
a steep red power spectrum. These induced residuals are 
correlated between different pulsar pairs as shown in Fig- 
ure [T] Although limits on the amplitude of the GWB can 
be obtained from the residuals of a single pulsar (see, e.g., 
iKaspiet al]|l994h . the GWB can only be detected with 
confidence by observing this pair-wise correlation. In pul- 
sar timing residuals, "red" (i.e., low-frequency) power can 
come from a variety of other physical effects. These in- 
clude irregular spindown behaviour known as "timing noise" 
iHobbs et al.ll2010d : [shannon fc Cordesll2010l . and references 
therein), vari ation in the pul se dispersion in the interstel- 
lar medium (|You et al] 20071) or cal ibration and other in- 
strumental errors l|van Straten|[200^ . There are also some 
sources of noise which are correlated between pulsars, such 
as instabilities in Terrestrial Time (TT) and inaccu racies in 
the Solar-System ephemeris l|Champion et al1l2010l ). An in- 
stability in TT will affect all pulsars in the same way, induc- 
ing a correlated signal which is independent of the angular 
separation of the pulsars on the sky, leading to a positive 
offset in the correlation curve in Figure [1] An inaccuracy in 
the Solar-System ephemeris will induce residuals which are 
positively correlated for pairwise angular separations less 
than 90 degrees. Such a signal could be correlated with the 
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GWB signal shown in Figure [T] All of these low-frequency 
variations are difficult to predict and need to be accounted 
for when implementing an algorithm to detect the GWB. 

3.1 The Expected GWB Signal 

Throughout this paper we assume a power-law form for the 
characteristic strain, hc{f), of th e isotropic stochastic GWB. 
This power-law is given by (e.g.. IPhinne vl l200ll :[ Jenet et al.l 
I2OO6I ) 

hcif)=A{f/f,y,r (2) 

where / is the GW frequency, /lyr = (lyr)~^ and A is a di- 
mensionless quantity termed the "amplitude" of the GWB. 
The smallest upper bound on the amplitude of the GWB 
from the literature is A < 1.1 x 10"" ijjenet et al.lf200(? ). 
The power-law form of the GWB is cons i stent with most 
models to date (e.g., Ijaffe fc Backer] l2003l : IWvithe fc Loebl 
\200A ). The spectral exponent, a, can take a range of val- 
ues depending on the source of the GWB under investiga- 
tion (e.g., cosmic strings, small-orbit black- ho le binaries). 
Howe ver, all predicted backgrounds have a < (|Jenet et al.l 
l2006l . and references therein), which results in a steeply de- 
creasing power spectrum in the timing residuals. Several 
models of the expected GWB from an ensemble of supermas- 
sive black-hole binaries (SMBHBs) predict that the ampli- 
tude of the GWB will b e in the range 5 x 10"^' ^ < < 10~" 
dJaffe fc Backer! l2003l : IWvithe fc Loebl l2003l : ISesana et al] 
|2008| ) with a spectral exponent a = 

If the GWB strain is described by equation ([2]) , then the 
pow er spectrum of t he induced ToA pertur bations is (see, 
e.g., lDetweilej[l979l : I Jenet et al.ll2005l . |2006| ) 

The cross power spectrum between the induced ToA pertur- 
bations in pulsars i and j is 

X^Jif) = PAf)aO^,) (4) 

where (J" (Oij) is given in equation [1] 

3.2 Detecting the GWB signal 

We estimate Xij{f) for each pair of pulsars. As the spectrum 
of the GWB is very steep, only the lowest frequencies are 
of interest. Fortunately, the irregular sampling has less ef- 
fect on the lower frequencies than on the higher frequencies 
because the low frequencies are heavily oversampled. The 
observations of each pair of pulsars overlap over some time 
span Toveriap. For Npsi = 20 there are A''pairs ~ 190 pairs. For 
each pair we estimate the cross power spectrum at harmonics 
of / = 1/Tovcriap. If the sampling were uniform, these esti- 
mates would be uncorrelated. In practice we find that they 
are not uncorrelated and this reduces the sensitivity of our 

1 ISesana et al.l l|2008l ) proposed a more complicated frequency- 
dependence for hc{f) involving an extra term proportional to 
f~^, which causes significant deviation from equation JJ} for 
/ > 10~* Hz; current PTA projects are not yet sensitive enough to 
distinguish between the two forms. Until a detection of the GWB 
is made, PTAs are expected to focus on frequencies / <C 10~* Hz. 



detection algorithm. It is probable that the independence 
can be restored using the Ch olesky spectr a l estim ation pro- 
cedure recently discussed by IColes et al.l (|2010t l . However, 
this is beyond the scope of this work. 

For some pairs, Toveriap can be much smaller than 
the length of one or both time series. For our time se- 
ries, Toveriap ranges from just 0.8 yr for PSRs J0437— 4715 
and J1824-2452, to 14.1 yr for PSRs J0711-6830 and 
J1939-I-2134. The use of only the overlapping residuals 
causes a bias in the cross power spectral estimates, the 
causes of which are currently not known. We correct this 
bias by removing a quadratic function from the overlapping 
section of the two time series using a weighted least-squares 
(WLSQ) fit, as shown in Figure [T] This fit is in addition 
to the standard timing model fit which estimates the pulsar 
parameters. We estimate the cross power spectrum: 

^.i(/)=^4/)-^;(/)/7^overlap (5) 

where J^i denotes the Discrete Fourier Transform (DFT) of 
the timing residuals of pulsar i and * denotes complex con- 
jugation. We use the following standard definition of the 
one-sided DFT: 

JV-l 

Hh) = 2 Yl ^e-^-^'=", (6) 

where i — %/— T in this particular case, N is the number of 
timing residuals, r„ is the n-th residual and k is an integer 
between 1 and (A^ — 1)/2, rounded down. Note that the k = 
term corresponds to the mean of the time series, which is 
zero for pulsar timing residuals. Calculating the DFT is not 
trivial because of the uneven sampling and variable error 
bars. We calculated J-i{fk) for every pulsar using a WLSQ 
fit of a sine term plus a cosine term at each fk = fc/Tovoriap- 
This gives identical results to a weighted Lomb-Scarglc esti- 
mate of the spectrum (Scargle 1982: Z echmeister fc Kiirsterl 
l2009l ). The variance of each cross-power spectral estimate is 

ai.,(/) = (P.(,f)>(^'.(,f))/2 (7) 

where (...) indicates an expectation value and Pi{f) is the 
spectral estimate of the residuals of pulsar i at frequency 
/. In practice, we calculate these expectation values using 
a power-law fit to the lowest frequencies in the spectrum of 
each pulsar. This power-law fit gives a spectral model for 
low frequencies in this pulsar. 

We account for the effects of fitting the timing model 
to the observations using two Monte Carlo simulations. The 
first simulation estimates the power spectrum — before and 
after pulsar parameter fitting — of simulated white noise 
with the same sampling and ToA errors as the residuals of 
each pulsar. Dividing the post-fit power spectrum by the pre- 
fit power spectrum gives the effective "tran sfer function" of 
the full TEMP02 fitting procedure (see, e.g.. lBlandford et al.l 
[1984; Helhngs 1989), and this process is repeated 1000 times 
to find the average transfer function. We describe this as an 
effective transfer function because the TEMP02 fitting pro- 
cess does not act exactly as a filter. We correct the measured 
cross power spectrum for each pulsar pair at each frequency 
by dividing by the geometric mean of the transfer functions 
of the two pulsars at that frequency. This correction iscom- 
mon between our analysis and that o f IVerbiest et"al] ()2009t ). 
but this is the only pulsar parameter fitting correction that 
Verbiest et al. perform. 
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However, the transfer function can only correct the ef- 
fects of the timing model fit as it acts on white noise in 
the residuals because, although fitting the timing model is 
a linear operation, it is not a filter. When the residuals are 
affected by red noise, fitting the full timing model to the 
residuals reduces P{f = 1/Toba) by considerably more than 
the white noise transfer function, where Tobs is the time span 
of the residuals. This is easily shown by simulation. A second 
correction is therefore necessary to measure the effect of the 
full timing model fit on the non-white GWB contribution to 
the residuals. We simulate ~10000 realisations of the resid- 
uals and add a simulated GWB signal with A — 'i x 10"^''' 
and a = —2 /3 to all pulsars using the method described in 
iHobbs et al.l ((2009 ). This value of A was chosen because it 
gives the largest GWB signal which is still small compared 
with the noise, hence reducing the number of required sim- 
ulations. We further reduce the number of simulations by 
fixing every pulsar to be at the same position and distance, 
giving the maximum correlated GWB signal between pul- 
sars. We perform the full pulsar parameter fit using tempo2, 
estimate the cross power spectrum in each realisation and 
apply the transfer function correction described above. We 
divide the average corrected cross power spectrum of each 
pulsar pair by the theoretical level of the cross power spec- 
trum given in equation Q. This process defines a set of 
"calibration factors", 7ij(/fc). When forming subsequent es- 
timates of the cross power spectrum using equation ((SJ, we 
calibrate each estimate at the lowest three frequencies of the 
cross power spectrum by dividing the cross-power-spectral 
estimate at frequency for pulsars i and j by Jij{fk)- 

After performing both of these corrections, we estimate 
A^. For each frequency channel, fk, of the cross power spec- 



trum (measured in yr 



have (cf. Equations [3] and 2]) 



= 127rVf"'"Real[X,,(A)] 



(8) 



where A~j indicates the measurement of A obtained from 
pulsars i and j and Real [Xij (fk)] is the real part of the cross 
power spectrum. The variance of Aij(^ (Gij) is then propor- 
tional to the variance of Xij . 

T o compare directly with the technique of Ijenet et al.l 
l|2005h . we perform a weighted sum of the Aij(^ {9ij) esti- 
mates over cross-spectral frequency to obtain a single esti- 
mate of Aij^ (Sij) for each pulsar pair. 



(9) 



where both summations range from k 
Nspcc,ij is the number of cross-spectral frequencies for pul- 
sars i and j. This final estimate of J^jC,{9ij) is similar to 
the unnormalised covariance between the residuals of pul- 
sars i and j. We also use the observed scatter in estimates 
of A^jC^{9ij) obtained from simulated observations to esti- 
mate the uncertainty SA^jC, {6ij) for each pulsar pair. 

Having fully calibrated our technique using simulations, 
we estimate the squared amplitude of the GWB, A^, by 
forming an average of the AljC,{6ij) estimates weighted by 
the inverse variance of each estimate, fn practice this average 
is done by performing a WLSQ fit to find the amplitude 
(and its corresponding uncertainty) for which the quantity 
A^C^ best fits the observed values of AljC,{6ij). For ease of 
notation, we index over all possible pulsar pairs using m. 



where m is an index running from 1 to A'paira and we set 
i^m = C {^ij)- In this case, the expression for is 



= 



and its unweighted variance is 



EmV^ 



2 



EmC/'^' 



(10) 



(11) 



This initial estimate of the error assumes that each of 
the SA^Cidij) are well-estimated. If this is not true, then 
we need to augment the error on A^ by an extra term which 
describes the amount of scatter in the residuals. This cor- 
responds to accounting for a non-unity reduced-x^ of the 
WLSQ fit which determines A^. Thus our final estimate for 
the variance of A^ is 



Em ( [^mCrr 



'A2 



Em ^"'/"^A^ 



1) 



Ei/*^: 



(12) 



which is just the weighted estimate of the variance of A^. If 
A"^ is significantly larger than cr^2i then a detection of the 
GWB has been made. This algorithm has been implemented 
as a TEMP02 plugirQ. 



4 RESULTS 

From the Verbiest et al. (2009) observations we estimate 
the squared GWB amplitude to be A^ = -4.5 x 10"^", with 
an uncertainty a ^2 = 9.1 x 10"'^". Our result is consistent 
with the null hypothesis, that there is no GWB present. 
Although the estimate is negative and therefore would lead 
to an unphysical GWB, it is not improbable because the 
standard deviation is a factor of 2 larger than the magnitude 
of the mean. We simulated many realisations of the Verbiest 
et al. (2009) observations, including the uncertainty given by 
the ToA error bars and a random process consistent with the 
low-frequency spectrum of the residuals but no GWB signal. 
These simulations showed that our estimate is consistent 
with the null hypothesis with 76 per cent confidence. This 
result is shown as a histogram in Figure [S] At first, one 
might think that this histogram could be used to provide a 
95 per cent confidence upper bound on the GWB amplitude. 
However, as discussed further below, any limit thus obta ined 
would not take account of self-noise (|jenet et al.ll2005l ) due 
to the GWB-induced perturbations at the pulsar. 

In Figure |4l we plot the 15 estimates of A\jC^(6ij) with 
the smallest uncertainty. It is clear from this figure that the 
current noise levels are larger than 4.5 x 10"'^° and our result 
is consistent with the null hypothesis. One might infer from 



the dot-dashed curve for A 



The 



1 X 10" 



that such a large 



C codes for this plugin and oth- 
ers from this paper are available at 
http://www.atnf .csiro.au/research/pulsar/teinpo2/index.php?n= 
Main. Plugins 
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Figure 3. The histogram shows the distribution of for simu- 
lations of the Verbiest et al. (2008, 2009) residuals with no GWB 
present. The thin dotted hne shows the value of obtained from 
the observations. The estimates to the right of the dotted line in- 
clude 76 per cent of the simulation results. All physical GWBs 
have > 0. 
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Figure 4. The 15 most precise estimates of ^ijC('^jj) the 
Verbiest et al. (2008, 2009) observations (points with error bars), 
the best-fit value of = —4.5 X 10"'^'' X f (dashed curve) and 
the signal expected from a strong GWB with A^ = 1 X 10~^* 
(dot-dashed curve). 



GWB signal is ruled out by the observations. These obser- 
vations may indeed rule out such a GWB signal, but if 
were actually 1 x 10~^* the noise levels on each Aij({9ij), 
which provide the upper bound, would be much higher. As 
the noise levels come from the power spectrum of the resid- 
uals of each pulsar, obtaining an upper bound using the 
noise levels is equivalent to obtaining an upper bound di- 
rectly from the individual power spectra and ignoring the 
cross correlations. We will not pursue this bounding tech- 
nique further in this paper as we are concentrating on the 
subject of detection. 



5 DISCUSSION 

The results of applying this algorithm to the Verbiest et al. 
(2008, 2009) data axe disappointing in the sense that the 
sensitivity is considerably lower than that calculated in the 
appendix to Verbiest et al. (2009) . We believe the estimated 
errors to be correct because they are calibrated by simula- 
tion, so we ask the question: Why are the cross power spec- 
tra of the GWB lower than expected? To investigate this we 
have run a series of simulation^ with GWB signals of differ- 
ing a mplitudes injected into the observations ( Hob bs et al.l 
l2009h . The results are shown in Figure [S] The mean values 
of the derived A"^ are plotted as solid lines connecting error 
bars (which indicate the uncertainty in the mean) for two 
cases: (1) the algorithm including correction with the ^ij 
calibration factors (thick solid line); and (2) the algorithm 
with 7ij = 1 (thin solid line). These results show that our 
method returns a GWB amplitude estimate A^^^ such that, 
on average, A^ut = ^fn- Figure |6] shows that this GWB sig- 
nal is at the correct level on average in every pulsar pair. 
The difference between the thick solid line and the thin solid 
line in Figure [5] indicates that the GWB power is reduced 
by a factor of ~12 because of the pulsar parameter fitting. 

We can estimate the amount of GWB signal lost in es- 
timation of different timing parameters by calculating the 
weighted average calibration factor in the lowest frequency 
channel of each pulsar pair. While this will be at a different 
frequency for each pair, it nevertheless provides a straight- 
forward figure of merit for comparing the effect of fitting 
different timing model parameters. For the full TEMP02 
fit acting on the Verbiest et al. (2008, 2009) residuals, we 
find 7i7(/ = 1/Tovcriap) = 0.0790 ± 0.0002, which repre- 
sents an average loss of 0.0790"^ — 12.7 in GWB signal at 
f — 1 /Toveriap ■ This explains the large decrease in sensitivity 
of our method compare d to that presented in the appendix 
of IVerbiest et al.l l|2009l ). which did not fully account for the 
effect of pulsar parameter estimation on the GWB signal. 
In Table [3] we show the weighted average calibration fac- 
tor at / = l/Toveriap whcn fitting for different parameters 
in the timing model. The estimation of the pulsar position 
and parallax have little effect on 7ij(/ = 1/Tovcriap) since 
Toveriap is a fcw timcs greater than 1 yr for most of our pul- 
sar pairs, and so are not shown in Table |3] This table in- 
dicates that one can almost determine the complete effect 
of fitting on the GWB sensitivity by only including fits for 
the spin frequency, its derivative and the arbitrary phase off- 
sets between different observing systems. Additionally, while 
the spin frequency derivative fit only significantly affects the 
power in the lowest frequency channel, the arbitrary phase 
offsets affect the power in the lowest few channels which can 
significantly affect our estimate of A^. 

The dashed lines in Figure [S] show that for GWB am- 
plitudes around ^4^ = 5 x 10""^", the average uncertainty on 



^ These simulations use a spread of pulsar distances and synthe- 
sise residuals with the same sampling as the Verbiest et al. (2008, 
2009) observations. The simulated residuals include white noise 
consistent with the observed error bars, red noise consistent with 
the spectral model mentioned in equation l[7) and a signal from 
a GWB with a = —2/3 and with a range of amplitudes between 
= 6.4 X 10~^^ and = 4 X 10~^*. We did not perform 
post-Kcplerian parameter fits. 
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Figure 5. Average as a function of input GWB for the 
Verbiest et al. (2008, 2009) residuals. The ordinate gives the av- 
erage output from our detection algorithm. The triple-dot- 
dashed line indicates points where the input is equal to the 
output A^. We have considered 2 cases: performing the full detec- 
tion procedure (thick lines) and the uncalibrated detection pro- 
cedure which uses ■yij(f) = 1 (thin lines). In both cases we have 
averaged over 1400 realisations for each input A^, and estimated 
the average output A^ (solid lines), where the error bars give the 
error in the mean of A^. The dashed lines give the square root 
of the average of cr^. in each case, and arc in good agreement 
with the sample standard deviations over the amplitude range of 
interest (dotted lines). 



is double the average uncertainty when there is no in- 
put GWB. This extra contribution to the uncertainty comes 
from the effect of the GWs passing near the pulsar, which 
we refer to as the "self-noise" of the GWB. For larger val- 
ues of A'^, the uncertaint y on is dominat ed by the GWB 
self-noise as discussed in IJenet et"aLl l|2005h . This provides 
a limitation on the confidence with which we can place an 
upper bound on the amplitude of the GWB. Because of the 
self-noise of the GWB, we can obtain at best an 80 per 
cent confidence upper bound on the GWB amplitude; we 
can never obtain a 95 per cent confidence bound with our 
current time series and weighting scheme. Furthermore, any 
limit obtained thus would be considerably worse than one 
obtained through other methods, such as direct power esti- 
mation, because of the huge variation in noise levels amongst 
our pulsarfl 

We confirm the accuracy of the measured uncertainty 
on each estimate of Aij( {6ij) using the reduced-x^ of the 
WLSQ fit that determines A'^. The reduced-x^ of this fit is 



(-^ pairs 



'Ala] - A:Xk 



(13) 



which has a value of 1.3 for the Verbiest et al. (2008, 2009) 



■* The lJenet et ah limit method requires that the timing 

residuals of each pulsar be white, so it cannot be used on these ob - 
servations. The method presented in Ivan Haasteren et al ] l|2009l ') 
could be applied to these observations, but this would require a 
large amount of computation time and any limit obtained would 
be difficult to confirm via Monte Carlo simulation. 



Figure 6. The expected covariance in simulated residuals which 
include a GWB component with squared amplitude A^ = 1 X 
'pjjg smooth dashed curve corresponds to the theoretical 
covariance for an input A^ = 1 X 10^^*. The points correspond 
to the mean estimates of ^fjC (^ij) (^^e Equations [l] and [2]| from 
200 simulated sets of timing residuals for the 20 PPTA pulsars. 
The error bars give the uncertainties in these mean estimates. For 
clarity we only plot the 20 pairs with the smallest rms scatter in 
their estimates of Af^C, (dij) over the 200 simulations. 



Timing 


Weighted 


Uncertainty 


Sensitivity 


Model 


mean of 


in Weighted 


Loss 


Parameters 


lijU = l/Toverlap) 


Mean 


Factor 


U, V 


0.1716 


0.0003 


5.83 


V, V, JUMP 


0.0796 


0.0002 


12.6 


ALL 


0.0790 


0.0002 


12.7 



Table 3. The effect of fitting different combinations of timing 
model parameters on the GWB signal in the lowest frequency 
channel. Values in the 4th column are the inverse of values in the 
2nd column. The symbols are: v (pulse frequency); v (pulse fre- 
quency derivative); "JUMP" (arbitrary phase offsets between dif- 
ferent observing systems were removed from all pulsars); "ALL" 
(all timing model parameters were fit). 



residuals, indicating that the uncertainty estimates cr^2 are 

consistent with the rms variation of the estimates A^. We 
obtain an independent estimate of the accuracy of the mea- 
sured errors by making use of the information contained 
in the imaginary part of the cross power spectrum, which 
we denote Imag [Xij(/)]. We calculate Imag [yl^jC^ (Sij)] 
by evaluating equation ((8| with Imag [Xij (/)] in place of 
Real [Xij(/)]. We then process Imag (^ij)] in exactly 

the same way as the real part is processed. Since correlation 
coefficients are real, we expect that Imag [AfjC^ (0ij)] will 
contain no correlated signal. This means that we can calcu- 
late the analogue of the reduced-x^ using Imag [^ijC (^ij)] • 



Xr.im 



1) 



E 



(Imag [Am 



(14) 



Similar to the reduced-x'^, if the errors on A^jC (dij) are well- 
estimated then this quantity should be near unity. For the 
Verbiest et al. (2008, 2009) residuals, we find x?,im = 0.87, 
indicating that the errors are well-estimated. 
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Table 4. The results from estimating with different estima- 
tors averaged over 10^ simulations of realistic residuals including 
a GWB with ^2 = 1 X 10~30. 



Estimator 


Mean 
(xlO-^O) 


Error in Mean 
A2 (x 10-30) 


rms of 
(xlO-^O) 


WLSQ [our method] 


0.99 


0.038 


12 


Unweighted LAD 


1.0 


0.038 


12 


Weighted LAD 


0.84 


0.041 


13 



Although both Xv ^nd Xr.im show that the uncertain- 
ties a^2 are reliable on average, these uncertainties come 
from power spectral estimates so they are random variables. 
We estimated the sensitivity of A'^ to variations in cr^2 by 

k 

multiplying each ct^2 by a random factor, distributed as the 

k 

square root of the product of two random variables with 
two degrees of freedom. This is the expected distribution for 
each o"a2- We found that a ^2 increased by a factor of 1.6, 
indicating that the use of incorrect SAij(^ (dij) estimates de- 
grades the sensitivity of the A'^ measurement by only a factor 
of 1.6. 

However, the ^i_,C(^»j) ^-re not Gaussian; rather they 
come from the sum of two pairwise products of independent 
Gaussian variables and thus have a two-sided exponential 
distribution which is reflected in Figure [3] This means that 
the maximum likelihood estimator for A^ is not a WLSQ 
estimator b ut a weigh ted least absolute deviation (LAD) 
fit (see, e.g.. [Cro3[200g V We tested both weighted and un- 
weighted LAD fits and found that the results for WLSQ and 
unweighted LAD fits were very similar, while the weighted 
LAD fit introduced a small bias in the mean. These results 
are shown in Table U We suspect that the bias occurs be- 
cause any LAD fit includes a 'dead-zone' feature, where a 
range of parameter estimates give the same minimum abso- 
lute deviation. This dead zone is negligible when the number 
of estimates is large, but can be significant otherwise. Since 
our A^ estimates are dominated by a small number of Aj. 
measurements and the results of the different estimators are 
similar, we chose the more standard WLSQ fit in calculating 
A'^. Although the WLSQ estimator is not maximum likeli- 
hood, it is apparently more robust in our particular case. 

Estimation of A^ is also largely independent of changes 
to the method of spectral analysis. We experimented with 
reducing the white noise in the residuals by smoothing each 
time series over a 60-day period before commencing the 
spectral analysis. We also tested interpolation using a con- 
strained cubic spline of each smoothed time series onto a 
14-day grid common to all pulsars before the spectral anal- 
ysis. The results of these different approaches are given in 
TableO Since there was no statistically significant difference 
between the different approaches, for simplicity we elected 
not to smooth or interpolate the residuals. 



5.1 Treatment of large GWB signals 

For their detection statistic. I Jenet et al.l (|2005l ) calculate the 
normalised cross correlation between the timing residuals 
of each pulsar pair. They optimise the S/N ratio using a 
filter designed to whiten the residuals before correlation. 



Table 5. The results from the Verbiest et al. (2008, 2009) obser- 
vations using different methods of spectral analysis of the timing 
residuals. 



Processing 


i2 


(7^2 


Performed 


(xlO-30) 


(xlO-30) 


Smoothing &i Interpolation 


3.0 


10 


Smoothing only 


-7.8 


10 


No smoothing [our method] 


-4.5 


9.1 



For a simulation of the 20 PPTA pulsars, this approach in- 
creased the maximum achievable detection significance for 
a GWB from 3cr to 13a. However, their filter cannot be 
applied to real pulsar timing observations without modifica- 
tion. We investigated the effect of such a filter by perform- 
ing simulations of the Verbiest et al. (2008, 2009) residuals 
where each simulation included a signal from a GWB with 
A > 3 X 10~^^. In the frequency domain, the filter takes the 
form of a weighting factor, so we optimised this weighting 
factor to match the large input GWB amplitude. We found 
that this method did not improve the S/N ratio, and we 
traced this under-performance to the problem of spectral 
leakage from the lowest frequencies to the higher frequen- 
cies. We found that the first few cross-spectral estimates, 
which make the largest contribution to our detection statis- 
tic, were all more than 90 per cent correlated with the lowest 
spectral estimate (i.e., at frequency / — l/Tovcriap), mean- 
ing that re- weighting cannot change the overall S/N ratio. 
The spectral leakage is particularly significant because of the 
irregular sampling and variable ToA uncertainties in these 
observations. We expect that an improved spectral analysis 
technique (e.g.. lColes et al.]|2010l ) will eliminate the spectral 
leakage and enable us to take advantage of more degrees of 
freedorr0 when the GWB signal is larger than the noise. 

5.2 Fitting timing models over different data 
spans 

The time series we consider in this paper have widely vary- 
ing time spans, which has not been a feature of most PTA 
analyses to date. As part of the pulsar parameter estimation, 
we fit for the pulse period and its derivative over the full du- 
ration of each time series. Originally, we then computed the 
cross power spectra from the overlapping portion of residu- 
als of each pulsar pair with no further processing. However, 
upon simulating this procedure, we found that the lowest 
frequencies in the cross power spectra were biased whenever 
Tobs > Tovcriap. This bias took the form of a significantly 
non-zero imaginary part in the cross power spectrum. Also, 
we found that much of the correlated signal at low frequen- 
cies was removed, as shown in Figure [T] We were unable to 
eliminate these effects unless we performed a WLSQ fit of 
a quadratic function for each time series over the overlap- 
ping time range. This restores the correlation in the GWB 

^ In contrast to lVerbiest et al"] ||2009|) which states that quadratic 
fitting removes one degree of freedom from the power spectrum 
of each pulsar's residuals, our analysis has shown that quadratic 
fitting does not affect the number of degrees of freedom in the 
lowest few frequency channels of each power spectrum. 
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Figure 7. The effect of fitting a timing model over different data 
spans. Tfie time series in tlie upper three panels are 5 years long, 
the time series in the lower three panels are 15 years long (to keep 
the y-axis scaling consistent, the plotting window has truncated 
the longer time series in the first 2 panels, and the bottom right 
panel only includes the overlapping data). The vertical dotted 
lines indicate the overlapping timing residuals for these time se- 
ries. We added the same large signal to both time series and the 
time series are identical in the overlapping region (left panels). 
After fitting the timing model (middle panels), this signal is no 
longer correlated between the two time series. The correlation is 
restored by performing a WLSQ fit of a quadratic function in the 
overlapping region of the two time series (right panels). 



independent of angular separation. Any estimate of the clock 
error will thus be correlated with the estimate of the GWB 
amplitude. Had we made a significant detection of the GWB, 
this would have to be accounted for. To estimate the impor- 
tance of possible clock instabilities, we processed the Ver- 
biest et al. (2008, 2009) observations us ing the ver sion of 
TT released by BIPM in 2010 (see, e.g., IPetitI 120031 ). This 
post-corrected timescale has revealed statistically significant 
inaccuracies in TT(TAI) . The results are shown in Table [51 
While the change of clock reference only changes our esti- 
mated GWB level by nine per cent of the uncertainty, the 
absolute change (0.8 x 10 ~^°) is at a significant level for some 
predictions of the GWB (jjaffe fc Backerll2003l : [ Sesana et al.l 
.2008) . This implies that such instabilities in TT must be 
accounted for when analysing future data sets. 

The results from using the newest Solar-System 
ephemeris DE421 (jPolkner et al.ll2009^ are given in Table |6l 
While there have been some improvements in this ephemeris 
version compared to DE405, most of the changes are ab- 
sorbed by the pulsar parameter fit. The estimated GWB 
level has changed by 24 per cent of the uncertainty. If we 
assume DE421 is correct, then the use of DE405 is similar 
to introducing a spurious GWB signal with A = 1.5 x 10~^^, 
a signal which is undetectable in most time series from the 
Verbiest et al. (2008, 2009) observations. However, future 
observations will need to account for the effects of inaccura- 
cies in the Solar-System ephemeris. 



Table 6. The results from using updated realisations of TT and 
the Solar-System ephemeris. The last column gives the change in 
the value of with respect to processing the observations with 
TT(TAI) and DE405, the realisations used for the Verbiest et al. 
(2008, 2009) residuals. 



Realisation 


Solar 






Change 


of Terrestrial 


System 






in 


Time 


Ephemeris 


(xlO-30) 


(xlO-30) 


(xlO-30) 


TT(TAI) 


DE405 


-4.5 


9.1 


0.0 


TT(TAI) 


DE421 


-2.3 


9.4 


2.2 


TT(BIPM2010) 


DE405 


-3.7 


8.7 


0.8 



signal between different pulsars (right panels of Figure [7|). 
This additional WLSQ fit will introduce a new bias because 
of removing some of the GWB signal at f — l/Toveriap, but 
this new bias is easily corrected with the calibration fac- 
tors 7ij(/). However, there is an additional loss of 10 per 
cent of the GWB signal in the Verbiest et al. (2008, 2009) 
observations because of this extra WLSQ fit. 

5.3 Correlated signals in the timing residuals 

The GWB analysis is complicated by the unknown effects 
of other correlated signals in the timing residuals. Instabil- 
ities in TT and errors in the Solar-System ephemeris both 
produce signals which are correlated between different pul- 
sars. We estimated the effect of these uncertainties by us- 
ing an updated timescale and the most recent Solar-System 
ephemeris. 

Instabilities in TT produce a positive cross correlation 



5.4 Contribution of different pulsars to A'^ 

It is difficult to determine the exact contributions to the 
weighting of each pulsar pair when using error bars derived 
from Monte Carlo simulations. The dominant effect is the 
size of Tovcriap • For a GWB caused by SMBHBs, the weight- 
ing factor increases approximately as Tovcriap ■ ^ higher noise 
level in the residuals of each pulsar in the pair will decrease 
the weight of that pair approximately linearly. The angle 
subtended at the observer by the pair of pulsars 9ij can be 
important if 9ij is near the zeroes of the function plotted in 
Figure [T] 

To determine which pulsars contribute the most to our 
estimate of the GWB, we perform the WLSQ fit described 
by Equations (fTU)) and (fTTj) to only 189 of the possible 190 
^ijC i^ij) estimates. By varying which estimate of j4^_,(" (Oij) 
is removed, we can find the pulsar pairs which have the 
greatest infiuence over the measurement of A'^ in these resid- 
uals. This is performed by finding AA'^ for each pair of pul- 
sars, which is the measured A"^ from all pulsar pairs minus 
the value of A'^ when not including the given pulsar pair. 
Those pairs with the largest contribution to this measure 
are given in Table [T] and a histogram of the absolute value 
|Aj42| for all pulsar pairs is provided in Figure |5] 

This analysis shows that the measurement of A'^ is de- 
termined by only a few pulsar pairs. This severely reduces 
the number of degrees of freedom when detecting the GWB, 
and thus decreases the maximum attainable detection con- 
fidence (see Jenet e t al. .2005i ) because it reduces our abil- 
ity to average out the self-noise in the residuals caused by 
the GWB signal at each pulsar. Observing more strong pul- 
sars is essential to increasing the number of degrees of free- 
dom in order to detect the GWB with reasonable confidence. 
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Table 7. The nine pulsar pairs whose absence from the array 
changes the measurement of from the Verbiest et al. (2008, 
2009) residuals by more than 1 X lO""^". The first column contains 
the names of the pulsars in the pair, the second column lists values 
of AA^, and the third column gives the change as a percentage 
of the value of derived when using all our data. 



Removed Pulsar Pair 


Ai2 (xlO-^O) 


Percentage change 


J1713+0747, J1744-1134 


18.0 


-400% 


J2124-3358, J2145-0750 


2.32 


-52% 


J1730-2304, J1744-1134 


2.10 


-47% 


J0711-6830, J2145-0750 


1.26 


-28% 


J0437-4715, J1909-3744 


-1.07 


24% 


J0437-4715, J2129-5721 


-1.36 


30% 


J0437-4715, J2145-0750 


-1.41 


31% 


J1713+0747, J2145-0750 


-3.97 


88% 


J0437-4715, J1713+0747 


-7.15 


159% 



10 10 10 
lAA^I estimate 



Figure 8. The effect on of the removal of different pulsar 
pairs, as measured by |Ay42|. Almost all pulsar pairs have no 
significant effect on the value of obtained from the Verbiest et 
al. (2008, 2009) residuals. 



This is further endo rsement of the Intern ational Pulsar Tim- 
ing Array concept ijHobbs et al.l[2010al ). but contrary to a 
suRgested st rategy for detection of individual GW sources 
iBurt et alllioiOl This is a fundamental difference between 
the single GW source detection problem and the GWB de- 
tection problem. 



6 CONCLUSIONS 

In implementing a GWB detection a l gorith m along the 
lines originally proposed bv ljenet et al.l (|2005l l we have con- 
fronted a number of issues which must be addressed when 
using real observations. We find that in practice the S/N 
ratio can be reduced by a fact or of ^12 com p ared with 
the ideal situation discussed by IVerbiest et all (j2009l ) be- 
cause of the fitting of a timing model to form the residuals. 
In particular, almost all of the signal loss is caused by the 
fitting of a quadratic term and arbitrary phase offsets be- 
tween different observing systems. We also find that it will 
be important to estimate and correct both clock errors and 
ephemeris errors when attempting to detect th e GWB at a 
level less than A — 2x 10~^^. As pointed out bv ljenet et al.l 



(|2005l ). prewhitening will be required to obtain detection sig- 
nificance larger than Scr. We find that this cannot be done 
without solving the problem of spectral leakage due to ir- 
regular sampling and variable ToA uncertainties. 

Fortunately, there are encouraging indications that 
many of these problems can be solve d. Recent work 
(iHobbs et al.ll2010dl : [Champion et~alll201(]| ') shows that clock 
errors and ephemeris errors can be estimated and removed. 
These errors are at a level which would disrupt the GWB 
signal in pulsar timing observations in the near future, and 
could even impact the analysis of a modified version of the 
Verbiest et al. (2008, 2009) observations which did not in- 
clude arbitrary phase offsets between observing systems. As 
systems with more sensitivity become available, the clock 
and ephemeris communities will improve their data sets. It 
appears possible to improve the process of fitting a timing 
model and also to im prove the sp e ctral leakage using the al- 
gorithm discussed bv lColes et al.l l|2010l) . It has proved pos- 
sible to calibrate most of the phase discontinuities between 
different observing systems in the PPTA observations and 
this alone can improve the S/N ratio by a factor of two. 

We have not discussed DM variations, but it is likely 
that some of the low frequency noise in our residuals is due to 
such interstellar propagation effects. Certainly as the various 
PTA data sets improve it will be essential to estimate and 
remove any frequency-dependent effects. 

Our analysis shows that, although the Verbiest et al. 
(2008, 2009) data set contains observations of 20 pulsars 
spanning many years, only a few of the pulsars in this data 
set contribute significantly to detecting the GWB, thereby 
reducing our detection confidence. It is uncertain whether 
this will be the case for the most recent observations from 
the PPTA. Observations of a larger sample of pulsars with 
precise ToA measurements will help to overcome this prob- 
lem. 
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